﻿// Copyright 2005, 2006 - Morten Nielsen (www.iter.dk)
//
// This file is part of SharpMap.
// SharpMap is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// SharpMap is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License
// along with SharpMap; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 

// SOURCECODE IS MODIFIED FROM ANOTHER WORK AND IS ORIGINALLY BASED ON GeoTools.NET:
/*
 *  Copyright (C) 2002 Urban Science Applications, Inc. 
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */
using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;
using Petra.Gis.ReferenceSystemServices;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{
  /// <summary>
  /// represents Universal (UTM) and Modified (MTM) Transverses Mercator projections.
  /// </summary>
  /// <remarks>
  /// <para>
  /// This is a cylindrical projection, in which the cylinder has been rotated 90°.
  /// Instead of being tangent to the equator (or to an other standard latitude),
  /// it is tangent to a central meridian. Deformation are more important as we
  /// are going futher from the central meridian. The Transverse Mercator
  /// projection is appropriate for region wich have a greater extent north-south
  /// than east-west.</para>
  /// 
  /// <para>Reference: John P. Snyder (Map Projections - A Working Manual,
  ///            U.S. Geological Survey Professional Paper 1395, 1987)</para>
  /// </remarks>
  [DataContract]
  public class TransverseMercatorProjection : ProjectionTransform
  {

    /* Variables common to all subroutines in this code file
  -----------------------------------------------------*/
    private double _semiMajor;		   //major axis
    //private double _semiMinor;		   //minor axis
    private double scale_factor;	   //scale factor
    private double central_meridian; //Center longitude (projection center)
    private double lat_origin;	     //center latitude
    private double e0, e1, e2, e3;	 //eccentricity constants
    private double e, es, esp;		   //eccentricity constants
    private double ml0;		           //small value m
    private double false_northing;	 //y offset in meters
    private double false_easting;	   //x offset in meters
    //static double ind;		/* spherical flag			*/




    /// <summary>
    /// xaml constr.
    /// </summary>
    public TransverseMercatorProjection() 
    { }





    /// <summary>
    /// Creates an instance of an TransverseMercatorProjection projection object.
    /// </summary>
    /// <param name="proj">ProjectedCoordinateSystem that must contain all necessary parameters</param>
    public TransverseMercatorProjection(ProjectedCoordinateSystem proj)
      : this(proj, false)
    {

    }





    /// <summary>
    /// Creates an instance of an TransverseMercatorProjection projection object.
    /// </summary>
    /// <param name="parameters">List of parameters to initialize the projection.</param>
    /// <param name="inverse">Flag indicating wether is a forward/projection (false) or an inverse projection (true).</param>
    /// <remarks>
    /// <list type="bullet">
    /// <listheader><term>Items</term><description>Descriptions</description></listheader>
    /// <item><term>semi_major</term><description>Semi major radius</description></item>
    /// <item><term>semi_minor</term><description>Semi minor radius</description></item>
    /// <item><term>scale_factor</term><description></description></item>
    /// <item><term>central meridian</term><description></description></item>
    /// <item><term>latitude_origin</term><description></description></item>
    /// <item><term>false_easting</term><description></description></item>
    /// <item><term>false_northing</term><description></description></item>
    /// </list>
    /// </remarks>
    public TransverseMercatorProjection(ProjectedCoordinateSystem proj, bool inverse)
      : base(proj, inverse)
    {
      //this.Name = "Transverse_Mercator";
      //this.Authority = "EPSG";
      //this.AuthorityCode = 9807;

      //ProjectionParameter par_semi_major = GetParameter("semi_major");
      //ProjectionParameter par_semi_minor = GetParameter("semi_minor");
      //ProjectionParameter par_scale_factor = GetParameter("scale_factor");
      //ProjectionParameter par_central_meridian = GetParameter("central_meridian");
      //ProjectionParameter par_latitude_of_origin = GetParameter("latitude_of_origin");
      //ProjectionParameter par_false_easting = GetParameter("false_easting");
      //ProjectionParameter par_false_northing = GetParameter("false_northing");

      ////Check for missing parameters
      //if (par_semi_major == null)
      //  throw new ArgumentException("Missing projection parameter 'semi_major'");
      //if (par_semi_minor == null)
      //  throw new ArgumentException("Missing projection parameter 'semi_minor'");
      //if (par_scale_factor == null)
      //  throw new ArgumentException("Missing projection parameter 'scale_factor'");
      //if (par_central_meridian == null)
      //  throw new ArgumentException("Missing projection parameter 'central_meridian'");
      //if (par_latitude_of_origin == null)
      //  throw new ArgumentException("Missing projection parameter 'latitude_of_origin'");
      //if (par_false_easting == null)
      //  throw new ArgumentException("Missing projection parameter 'false_easting'");
      //if (par_false_northing == null)
      //  throw new ArgumentException("Missing projection parameter 'false_northing'");

    }





    ///// <summary>
    ///// TransverseMercator projection axis list
    ///// </summary>
    //public override AxisList Axes
    //{
    //  get
    //  {
    //    // this is default TransverseMercator axis list with correct orientation
    //    // this orientation must remain as defined because of projection equations
    //    return new AxisList()
    //    {
    //      new Axis("X", Axis.East),
    //      new Axis("Y", Axis.North)
    //    };
    //  }
    //}





    protected override bool parametersValid()
    {
      return Parameters.ContainsKey("scale_factor")
          && Parameters.ContainsKey("central_meridian")
          && Parameters.ContainsKey("latitude_of_origin")
          && Parameters.ContainsKey("false_easting")
          && Parameters.ContainsKey("false_northing");

    }


    



    protected override void updateTransform()
    {
      //if (!(_projCS.Unit is DistanceUnit))
      //  throw new ArgumentException(this.ToString() + " projection must have DistanceUnit");

      //_projectedUnit = (DistanceUnit)_projCS.Unit;

      scale_factor = Parameters["scale_factor"];
      central_meridian = Parameters["central_meridian"] * AngularParameterUnit.RadiansPerUnit;
      lat_origin = Parameters["latitude_of_origin"] * AngularParameterUnit.RadiansPerUnit;
      false_easting = Parameters["false_easting"] * DistanceParameterUnit.MetersPerUnit;
      false_northing = Parameters["false_northing"] * DistanceParameterUnit.MetersPerUnit;

      _semiMajor = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMajorAxis
        * Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      //_semiMinor = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMinorAxis
      //  * Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      //es = 1.0 - Math.Pow(_semiMinor / _semiMajor, 2);
      //e = Math.Sqrt(es);
      e = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.Eccentricity;
      es = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.e2;

      e0 = ProjMath.E0fn(es);
      e1 = ProjMath.E1fn(es);
      e2 = ProjMath.E2fn(es);
      e3 = ProjMath.E3fn(es);
      ml0 = _semiMajor * ProjMath.Mlfn(e0, e1, e2, e3, lat_origin);
      esp = es / (1.0 - es);
    }





    /// <summary>
    /// Converts coordinates in decimal degrees to projected meters.
    /// </summary>
    /// <param name="lonlat">The point in radians</param>
    /// <returns>Point in meters</returns>
    protected override void geographicToProjectedN(ICoordinate lonlat)
    {

      //// convert to radians
      //lonlat.lon *= _geogUnit.RadiansPerUnit;
      //lonlat.lat *= _geogUnit.RadiansPerUnit;

      double delta_lon = 0.0;	 // Delta longitude (Given longitude - center
      double sin_phi, cos_phi; // sin and cos value
      double al, als;		       // temporary values
      double c, t, tq;	       // temporary values
      double con, n, ml;	     // cone constant, small m

      delta_lon = ProjMath.AdjustLongitude(lonlat.lon - central_meridian);
      ProjMath.Sincos(lonlat.lat, out sin_phi, out cos_phi);

      al = cos_phi * delta_lon;
      als = Math.Pow(al, 2);
      c = esp * Math.Pow(cos_phi, 2);
      tq = Math.Tan(lonlat.lat);
      t = Math.Pow(tq, 2);
      con = 1.0 - es * Math.Pow(sin_phi, 2);
      n = _semiMajor / Math.Sqrt(con);
      ml = _semiMajor * ProjMath.Mlfn(e0, e1, e2, e3, lonlat.lat);
      
      lonlat.X = scale_factor * n * al * (1.0 + als / 6.0 * (1.0 - t + c + als / 20.0 *
        (5.0 - 18.0 * t + Math.Pow(t, 2) + 72.0 * c - 58.0 * esp))) + false_easting;
      
      lonlat.Y = scale_factor * (ml - ml0 + n * tq * (als * (0.5 + als / 24.0 *
        (5.0 - t + 9.0 * c + 4.0 * Math.Pow(c, 2) + als / 30.0 * (61.0 - 58.0 * t
        + Math.Pow(t, 2) + 600.0 * c - 330.0 * esp))))) + false_northing;

      //// convert to projcs units
      //lonlat.X /= _projUnit.MetersPerUnit;
      //lonlat.Y /= _projUnit.MetersPerUnit;
    }





    /// <summary>
    /// Converts coordinates in projected meters to decimal degrees.
    /// </summary>
    /// <param name="p">Point in meters</param>
    /// <returns>Transformed point in radians</returns>
    protected override void projectedToGeographicN(ICoordinate proj)
    {
      double con, phi;		/* temporary angles				*/
      double delta_phi;	/* difference between longitudes		*/
      long i;			/* counter variable				*/
      double sin_phi, cos_phi, tan_phi;	/* sin cos and tangent values	*/
      double c, cs, t, ts, n, r, d, ds;	/* temporary variables		*/
      long max_iter = 6;			/* maximun number of iterations	*/

      //// convert from projcs to meters
      //proj.X *= _projUnit.MetersPerUnit;
      //proj.Y *= _projUnit.MetersPerUnit;

      //todo: update units of false_easting and false_northing here
      proj.X -= false_easting;
      proj.Y -= false_northing;

      con = (ml0 + proj.Y / scale_factor) / _semiMajor;
      phi = con;
      for (i = 0; ; i++)
      {
        delta_phi = ((con + e1 * Math.Sin(2.0 * phi) - e2 * Math.Sin(4.0 * phi) + e3 * Math.Sin(6.0 * phi))
          / e0) - phi;
        phi += delta_phi;
        if (Math.Abs(delta_phi) <= ProjMath.EPSLN) break;
        if (i >= max_iter)
          throw new ApplicationException("Latitude failed to converge");
      }

      if (Math.Abs(phi) < ProjMath.HALF_PI)
      {
        ProjMath.Sincos(phi, out sin_phi, out cos_phi);
        tan_phi = Math.Tan(phi);
        c = esp * Math.Pow(cos_phi, 2);
        cs = Math.Pow(c, 2);
        t = Math.Pow(tan_phi, 2);
        ts = Math.Pow(t, 2);
        con = 1.0 - es * Math.Pow(sin_phi, 2);
        n = _semiMajor / Math.Sqrt(con);
        r = n * (1.0 - es) / con;
        d = proj.X / (n * scale_factor);
        ds = Math.Pow(d, 2);

        proj.lat = phi - (n * tan_phi * ds / r) * (0.5 - ds / 24.0 * (5.0 + 3.0 * t +
          10.0 * c - 4.0 * cs - 9.0 * esp - ds / 30.0 * (61.0 + 90.0 * t +
          298.0 * c + 45.0 * ts - 252.0 * esp - 3.0 * cs)));
        proj.lon = ProjMath.AdjustLongitude(central_meridian + (d * (1.0 - ds / 6.0 * (1.0 + 2.0 * t +
          c - ds / 20.0 * (5.0 - 2.0 * c + 28.0 * t - 3.0 * cs + 8.0 * esp +
          24.0 * ts))) / cos_phi));
      }
      else
      {
        proj.lat = ProjMath.HALF_PI * ProjMath.Sign(proj.Y);
        proj.lon = central_meridian;
      }

      //// convert to geogcs units
      //proj.lon /= _geogUnit.RadiansPerUnit;
      //proj.lat /= _geogUnit.RadiansPerUnit;
    }





    /// <summary>
    /// Returns the inverse of this projection.
    /// </summary>
    /// <returns>IMathTransform that is the reverse of the current projection.</returns>
    public override SRSTransform Inverse()
    {
      return new TransverseMercatorProjection(Projection, !IsInverse);
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }





    public static string[] ProjectionNames
    {
      get { return new string[] { "transverse_mercator", "transverse mercator" }; }
    }





    //new public static TransverseMercatorProjection Create(ProjectedCoordinateSystem pcs)
    public static SRSTransform Create(ProjectedCoordinateSystem pcs, bool getAsInverse)
    {
      return new TransverseMercatorProjection(pcs, getAsInverse);
    }


  }
}
